/*****************************************************************************FILE:  qrsfilt.cppAUTHOR:	Patrick S. HamiltonREVISED:	5/13/2002  ___________________________________________________________________________qrsfilt.cpp filter functions to aid beat detecton in electrocardiograms.Copywrite (C) 2000 Patrick S. HamiltonThis file is free software; you can redistribute it and/or modify it underthe terms of the GNU Library General Public License as published by the Free
Software Foundation; either version 2 of the License, or (at your option) any
later version.

This software is distributed in the hope that it will be useful, but WITHOUT ANY
WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
PARTICULAR PURPOSE.  See the GNU Library General Public License for more
details.

You should have received a copy of the GNU Library General Public License along
with this library; if not, write to the Free Software Foundation, Inc., 59
Temple Place - Suite 330, Boston, MA 02111-1307, USA.

You may contact the author by e-mail (pat@eplimited.edu) or postal mail
(Patrick Hamilton, E.P. Limited, 35 Medford St., Suite 204 Somerville,
MA 02143 USA).  For updates to this software, please visit our website
(http://www.eplimited.com).
  __________________________________________________________________________
	This file includes QRSFilt() and associated filtering files used for QRS	detection.  Only QRSFilt() and deriv1() are called by the QRS detector	other functions can be hidden.	Revisions:		5/13: Filter implementations have been modified to allow simplified			modification for different sample rates.*******************************************************************************/#include <math.h>#include "qrsdet.h"// Local Prototypes.int lpfilt( int datum ,int init) ;int hpfilt( int datum, int init ) ;int deriv1( int x0, int init ) ;int deriv2( int x0, int init ) ;int mvwint(int datum, int init) ;/******************************************************************************* Syntax:*	int QRSFilter(int datum, int init) ;* Description:*	QRSFilter() takes samples of an ECG signal as input and returns a sample of*	a signal that is an estimate of the local energy in the QRS bandwidth.  In*	other words, the signal has a lump in it whenever a QRS complex, or QRS*	complex like artifact occurs.  The filters were originally designed for data*  sampled at 200 samples per second, but they work nearly as well at sample*	frequencies from 150 to 250 samples per second.**	The filter buffers and static variables are reset if a value other than*	0 is passed to QRSFilter through init.*******************************************************************************/int QRSFilter(int datum,int init)	{	int fdatum ;	if(init)		{		hpfilt( 0, 1 ) ;		// Initialize filters.		lpfilt( 0, 1 ) ;		mvwint( 0, 1 ) ;		deriv1( 0, 1 ) ;		deriv2( 0, 1 ) ;		}	fdatum = lpfilt( datum, 0 ) ;		// Low pass filter data.	fdatum = hpfilt( fdatum, 0 ) ;	// High pass filter data.	fdatum = deriv2( fdatum, 0 ) ;	// Take the derivative.	fdatum = abs(fdatum) ;				// Take the absolute value.	fdatum = mvwint( fdatum, 0 ) ;	// Average over an 80 ms window .	return(fdatum) ;	}/**************************************************************************  lpfilt() implements the digital filter represented by the difference*  equation:** 	y[n] = 2*y[n-1] - y[n-2] + x[n] - 2*x[t-24 ms] + x[t-48 ms]**	Note that the filter delay is (LPBUFFER_LGTH/2)-1***************************************************************************/int lpfilt( int datum ,int init)	{	static long y1 = 0, y2 = 0 ;	static int data[LPBUFFER_LGTH], ptr = 0 ;	long y0 ;	int output, halfPtr ;	if(init)		{		for(ptr = 0; ptr < LPBUFFER_LGTH; ++ptr)			data[ptr] = 0 ;		y1 = y2 = 0 ;		ptr = 0 ;		}	halfPtr = ptr-(LPBUFFER_LGTH/2) ;	// Use halfPtr to index	if(halfPtr < 0)							// to x[n-6].		halfPtr += LPBUFFER_LGTH ;	y0 = (y1 << 1) - y2 + datum - (data[halfPtr] << 1) + data[ptr] ;	y2 = y1;	y1 = y0;	output = y0 / ((LPBUFFER_LGTH*LPBUFFER_LGTH)/4);	data[ptr] = datum ;			// Stick most recent sample into	if(++ptr == LPBUFFER_LGTH)	// the circular buffer and update		ptr = 0 ;					// the buffer pointer.	return(output) ;	}/*******************************************************************************  hpfilt() implements the high pass filter represented by the following*  difference equation:**	y[n] = y[n-1] + x[n] - x[n-128 ms]*	z[n] = x[n-64 ms] - y[n] ;**  Filter delay is (HPBUFFER_LGTH-1)/2******************************************************************************/int hpfilt( int datum, int init )	{	static long y=0 ;	static int data[HPBUFFER_LGTH], ptr = 0 ;	int z, halfPtr ;	if(init)		{		for(ptr = 0; ptr < HPBUFFER_LGTH; ++ptr)			data[ptr] = 0 ;		ptr = 0 ;		y = 0 ;		}	y += datum - data[ptr];	halfPtr = ptr-(HPBUFFER_LGTH/2) ;	if(halfPtr < 0)		halfPtr += HPBUFFER_LGTH ;	z = data[halfPtr] - (y / HPBUFFER_LGTH);	data[ptr] = datum ;	if(++ptr == HPBUFFER_LGTH)		ptr = 0 ;	return( z );	}/******************************************************************************  deriv1 and deriv2 implement derivative approximations represented by*  the difference equation:**	y[n] = x[n] - x[n - 10ms]**  Filter delay is DERIV_LENGTH/2*****************************************************************************/int deriv1(int x, int init)	{	static int derBuff[DERIV_LENGTH], derI = 0 ;	int y ;	if(init != 0)		{		for(derI = 0; derI < DERIV_LENGTH; ++derI)			derBuff[derI] = 0 ;		derI = 0 ;		return(0) ;		}	y = x - derBuff[derI] ;	derBuff[derI] = x ;	if(++derI == DERIV_LENGTH)		derI = 0 ;	return(y) ;	}int deriv2(int x, int init)	{	static int derBuff[DERIV_LENGTH], derI = 0 ;	int y ;	if(init != 0)		{		for(derI = 0; derI < DERIV_LENGTH; ++derI)			derBuff[derI] = 0 ;		derI = 0 ;		return(0) ;		}	y = x - derBuff[derI] ;	derBuff[derI] = x ;	if(++derI == DERIV_LENGTH)		derI = 0 ;	return(y) ;	}/****************************************************************************** mvwint() implements a moving window integrator.  Actually, mvwint() averages* the signal values over the last WINDOW_WIDTH samples.*****************************************************************************/int mvwint(int datum, int init)	{	static long sum = 0 ;	static int data[WINDOW_WIDTH], ptr = 0 ;	int output;	if(init)		{		for(ptr = 0; ptr < WINDOW_WIDTH ; ++ptr)			data[ptr] = 0 ;		sum = 0 ;		ptr = 0 ;		}	sum += datum ;	sum -= data[ptr] ;	data[ptr] = datum ;	if(++ptr == WINDOW_WIDTH)		ptr = 0 ;	if((sum / WINDOW_WIDTH) > 32000)		output = 32000 ;	else		output = sum / WINDOW_WIDTH ;	return(output) ;	}